LIBRARY 

NAVAL   P©STGRADUATE   SCHOOL 
MONTEREY.  CALIFORNIA    9394Q 


NPS-61Gr73041A 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


6  April  1973 


NPS-6lGr73041A 


Bubble  Transport  Theory  With  Application 
to  the  Upper  Ocean 


G.A.  Garrettson 


Physics  &  Chemistry 
Department 


Approved  for  public  release;  distribution  unlimited, 


FEDDOCS 

D  208.14/2 

NPS-61GR73041A 


10?.  I4l^  , 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


DUDLEY  KNOX  LIBRARY 
NAVAL  POSTGRADUATE  SCH' 
MONTEREY  CA  93943-5101 


Rear  Admiral  M.  B.  Freeman  M.  U.  Clauser 

Superintendent  Provost 

TITLE:       Bubble  Transport  Theory  with  Application  to  the  Upper  Ocean  * 

AUTHOR:     Garrett  A.  Garrettson 

ABSTRACT: 

The  formalism  of  transport  theory  is  adapted  to  a  general  description  of 
bubble  populations  in  a  moving  fluid.  The  bubble  distribution,  as  a  function 
of  position,  velocity,  radius,  and  time,  satisfies  a  Boltzmanntype  transport 
equation  that  is  derived  and  then  formally  solved  by  the  method  of  character- 
istics. Before  this  new  analytical  tool  can  be  applied,  properties  of  the 
medium  must  be  specified  and  a  bubble  dynamics  model  must  be  chosen.  General 
expressions  are  written  for  bubble  acceleration  and  radius  change  rate,  and 
known  models  of  bubble  gas  diffusion  and  drag  are  summarized  for  gas  bubbles 
in  liquids.  Application  to  the  upper  ocean  is  discussed  and  then  illustrated 
with  seme  sample  calculations.  Explicit  solutions  are  written  for  a  steady- 
state,  one  dimensional  ocean  with  distributed  sources;     the  results  clarify 
relations  between  observed  bubble  populations,  proposed  bubble  source  mechanisms, 
and  known  models  of  single  bubble  dynamics. 

This  research  was  supported  by:  Naval  Ship  Systems  Command,  Code  (PMS  302)  . 

Approved  by: 


TABLE  OF  CONTENTS 

I .  INTRODUCTION 5 

II.  THE  BUBBLE  TRANSPORT  EQUATION 7 

III.  SOLUTION  BY  METHOD  OF  CHARACTERISTICS 13 

A.  LAGRANGIAN  AND  EULERIAN  REPRESENTATIONS 14 

B.  FORMAL  SOLUTION 15 

IV.  SINGLE  BUBBLE  DYNAMICS 19 

A.  ACCELERATION 19 

B.  DRAG  ON  GAS  BUBBLES  IN  A  LIQUID 20 

C.  RADIUS  CHANGE  RATE 21 

D.  GAS  DIFFUSION 22 

V.  APPLICATION  TO  THE  UPPER  OCEAN 25 

A.  BUBBLE  SOURCES  AND  FLUID  VELOCITY  FIELDS 25 

B.  EXAMPLE 27 

1.  Characteristic  Curves 29 

2.  Bubble  Distributions 34 

3.  Numerical  Results 38 

VI .  CONCLUSION 42 

APPENDIX 43 

REFERENCES 46 

FIGURE  CAPTIONS 47 

FIGURES 49 

INITIAL  DISTRIBUTION  LIST 61 

FORM  DD  1473 65 


ACKNOWLEDGEMENTS 

The  author  would  like  to  thank  several  people  who 
aided  this  research  effort:  Professor  Herman  Medwin 
supplied  valuable  suggestions  and  information  concerning 
bubble  phenomena  in  the  upper  ocean;  LT.  Frank  H.  Heistand, 
USN,  aided  in  the  preliminary  investigations  of  single 
bubble  dynamics  and  provided  programming  support; 
LT.  Thomas  C.  Vajda,  USN,  supplied  programming  support  for 
obtaining  the  numerical  results  of  Section  V;  Mrs.  Clarey 
and  Mrs.  Sutherland  (NPS)  diligently  typed  the  manuscript; 
and  my  wife  Gloria  typed  corrections  and  gave  the  necessary 
moral  support. 

This  research  was  partially  supported  by  the  Naval 
Ship  Systems  Command,  Code  (PMS  302). 


I .   INTRODUCTION 
The  solutions  of  many  experimental  and  engineering  problems  depend  upon  a 
knowledge  of  bubble  distributions  in  a  fluid  medium.  For  instance,  gas  bubbles 
near  the  surface  of  the  ocean  are  important  in  a  variety  of  subjects  that 
include  underwater  sound  propagation  (Shulkin  1968  &  1969) ,  meteorology 
(Blanchard  &  Woodcock. 1957) ,  sea  surface  chemistry  (Sutcliffe,  et  al  1963), 
cavitation  (Fox  and  Herzfeld  1954) ,  and  air-sea  gas  exchange  (Kanwisher 

1963) .  Several  efforts  have  been  made  to  measure  bubble  densities  in 
the  ocean  (Blanchard  &  Woodcock  1957;  Glotov,  et  al  1962;  McCartney  and  Bary 
1965) ,  the  most  recent  made  by  Medwin  (1970) .  To  infer  near  surface  bubble 
populations  from  in-situ  measurements  of  acoustic  attenuation,  Medwin  exploited 
the  fact  that  a  gas  bubble's  scattering  and  absorption  cross  sections  for  sound 
at  the  bubble  resonant  frequency  can  typically  be  1000  times  its  geometrical 
cross  section.  For  this  reason  large  bubble  densities  near  the  sea  surface  can 
significantly  affect  underwater  sound  transmission. 

Many  complex  factors  affect  bubble  distributions  in  a  fluid  such  as  the 
ocean,  and  some  of  these  phenomena  have  been  carefully  investigated.  For 
example,  LeBlond  (1969)  examined  gas  diffusion  from  an  ascending  bubble,  and 
Levich  (1962,  Ch.  8)  is  an  excellent  reference  on  single  bubble  dynamics.  The 
purpose  of  this  paper  is  to  develop  and  to  demonstrate  new  analytical  tools 
for  investigating  relationships  between  observed  bubble  distributions  and  the 
phenomena  which  create,  transport,  and  eliminate  bubbles. 

In  sections  two  and  three  the  formalism  of  transport  theory  is  adapted  to 
the  general  problem  of  describing  bubble  populations  in  a  moving  fluid.  The 
bubble  distribution,  as  a  function  of  position,  time,  velocity,  and  radius, 
satisfies  a  Boltzmann-type  transport  equation  that  is  derived.  Using  a 
Lagrangian  viewpoint,  the  formal  solution  is  expressed  as  a  line  integral  of 
the  volume  source  function  along  characteristic  curves.  In  section  four  general 


expressions  are  obtained  for  bubble  acceleration  and  radius  change  rate  and 
sane  known  models  of  gas  diffusion  and  bubble  drag  are  introduced.  Finally, 
in  the  context  of  gas  bubble  transport  in  the  upper  ocean,  a  simple  model  is 
chosen  to  illustrate  the  utility  of  transport  theory.  Calculated  results 
clarify  relationships  between  observed  bubble  distributions,  proposed  bubble 
source  mechanisms,  and  known  models  of  single  bubble  dynamics. 


II.   THE  BUBBLE  TRANSPORT  EQUATION 

Consider  the  problem  of  describing  the  bubble  distribution  —  as  a  function 
of  position,  r,  velocity ,  v,  size,  I,   and  time^t—  in  a  moving  fluid  that  may 
contain  bubble  sources  and  sinks.  In  general,  a  "bubble"  might  be  considered 
any  simple,  closed  region  containing  a  fluid  somehow  different  from  that  of 
the  transporting  medium.  Since  examples  chosen  in  this  paper  refer  primarily 
to  spherical  gas  bubbles,  I  will  represent  the  bubble  radius.  However,  if  the 
bubbles  are  not  spherical,  then  £  might  represent  the  radius  of  a  spherical 
volume  equal  to  the  actual  bubble  volume,  or  one  might  wish  to  generalize  the 
description  to  include  more  than  one  parameter  for  size. 

Suppose  that  the  time  evolution  of  the  position,  the  velocity, 
and  the  radius  of  each  bubble  depends  only  upon  the  properties 
of  the  medium  and  the  bubble  characteristics  and  is  statistically  independent 
of  the  state  of  any  other  bubble.  Then  define  the  bubble  distribution  function, 
iMr,v,£,t)d3rd3vd£,  as  the  mean  number  of  bubbles  at  time  t  in  the  volume  d3r 
about  r,  with  velocity  in  d^v  about  v  and  the  radius  in  d£  about  l.*     Further- 
more, if  bubbles  interact  with  one  another,  suppose  that  no  more  than  two 
interact  at  any  one  time  and  that  the  collision  time  is  short  compared  to  the 
time  for  \\>   to  change  appreciably.  Then  \\>   satisfies  a  Boltzmann-type  transport 
equation  whose  solution  describes  the  ensemble  average  behavior  of  the  bubble 
population  in  terms  of  single  bubble  behavior. 

At  this  point  it  is  convenient  to  define  the  seven-dimensional  volume  ele- 
ment d7x  =  d3rd3vd£  and  to  adopt  a  Lagrangian  viewpoint.  In  the  time  interval 
(t,t+dt)  the  bubbles  that  were  in  d7x  about  (r,v,£)  have  moved  to  a  new  volume 


*E.g.,  in  the  Cartesian  representation,  position  components  lie  in  the  intervals 
(x,x+dx) ,  (y,y+dy) ,  and  (z,z4dz) ;  velocity  components  lie  in  the  intervals 
(vx,vx+dvx) ,  (Vy,vy+dVy) ,  and  (vz,vz4dvz);  and  the  radius  is  in  the  interval 
(J£,JL4d£). 
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d7x  about  (r',v ',£'),  where  the  primed  quantities  refer  to  the  time  t4dt  and 
the  unprimed  quantities  refer  to  the  time  t.  This  transport  occurs  by  a  change 
in  position  due  to  a  velocity,  v  =  dr/dt,  by  a  change  in  velocity  due  to  an 
acceleration,  a  =  dv/dt,  and  by  a  change  in  radius  due  to  a  rate  of  change, 
v  e  dA/dt: 

r'  =  r  +  vdt  +  0(dt2) 

v1  =  v  +  adt  +  0(dt2)  (1) 

I*   =  i  +  vdt  +  0(dt2) 

(Taylor's  series  about  time  t.)  For  gas  bubbles  in  the  ocean  and  elsewhere, 
the  acceleration  results  from  several  influences  that  include  gravity,  buoyancy, 
and  entrainment  of  bubbles  by  the  fluid,  while  the  rate  of  change  of  bubble 
radii  is  caused  by  hydrodynamic  compression,  surface  tension  and  gas  diffusion. 
Expressions  for  a  and  v  are  obtained  in  section  four  where  single  bubble 
dynamics  is  discussed. 

Conservation  of  bubbles  demands  that  the  number  of  bubbles  in  the  volume 
element  at  time  tfdt,  ^(r1 ,v' ,£' ,t+dt)d7x' ,  equals  the  number  that  start  in  the 
volume  element  at  time  t,  ^(r,v, £,t)d7x,  plus  (minus)  any  bubbles  that  are 
introduced  (lost)  by  distributed  sources  (sinks),  S(r,v,£,t)d/xdt: 

^(r'  ,v'  ,l'  ,t+at)a7T'  -   iMr,v,£,t)d7T  =  S(r,v,£,t)d7xdt  +  0(dt2),   (2) 

where  S(r,v,  £t)d7x  is  the  rate  at  which  bubbles  are  introduced  into  d7x  about 
(r,v,£)  by  distributed  sources.  To  develop  equation  (2)  it  is  sufficient  to 
note  that  terms  of  0(dt2)  will  disappear  in  the  limit  dt  ■+■  o.  Therefore  we 
expand  ijj1  =  ^(r1  ,v'  ,V  ,t4dt)  in  a  Taylor  series  about  ip  e  ip(r,v,Jl,t)  and  keep 

8 


terms  through  0  (dt) : 


*'  =  *  +   (3i|*/3t)dt  +   (v-Vr^)dt  +   (a.V^)dt  +  v(3^/H)dt  +  0(dt2)        (3) 


where  7r(.)  e  (3  (-)/3x,3(-)/3y,3(«)/3z)    and  V    (.)    e    (3(-)/3v  ,3(«)/3v  ,3(-)/3v  ) 

are  Cartesian  representations  of  the  "del"  operators  in  position  space  and 
velocity  space,  respectively. 

To  proceed,  d7r'  must  be  expressed  in  terms  of  d7x: 

dV  =  y(t+dt/t)d7T  (4a) 

where 

3(r,,vl,£') 


Y(t-+dt,t)    i   dgjfi)  idet   ((  q\/  q±))  (4b) 

l/J   =   1,2,...,/ 


is  the  Jacobian  with 

q±   e  x,  q2  e  y,  q3  e  z,  q4  =  v^,  q5  e  v  ,  qg  =  v^   and  q?  =  l;  (5) 

4^'-=  x'/  q2'  -  y,/  etc# 

The  expression  derived  for  y   in  the  appendix  can  be  expanded  in  a  Taylor  series 
about  t,  or  equation.  (1)  can  be  used  to  compute  the  determinant  in  (4b)  directly. 
In  either  case  the  result  is 

y(t+dt,t)  =  1  +  (V  -a)dt  +  (3v/3S,)dt  +  0(dt2)  (6) 

Substituting  (3) ,  (4)  and  (6)  into  (2) ,  collecting  terms  of  0  (dt) ,  dividing  by 
d7xdt,  and  taking  the  limit  dt  -*  o,  yields  the  bubble  transport  equation: 


->  -> 


3^/at  +  v-v  ifi  +  a«v  ip  +  vdty/dl  =  S-z  i>  (7) 

where  E  =  V  -a  +  3v/3£  (8) 

and  where      -°°  <  x,y,z  <  <*>,   -°°  <  v  ,v  ,v  <  »,  and  o  <  5,  <  °°. 

A  simpler  (bat  less  rigorous)  derivation  of  equation  (7)  is  obtained  by 
recognizing  that  vtJj,  a^,  and  v^  are  the  fluxes  of  bubbles  in  position  space, 
velocity  space,  and  radius  space,  respectively.  Then 

F  e  -  {V  •  (v\|>)  +  V  .  (aip)  +  3(vip)/3£} 

is  the  net  flux  of  bubbles  into  d7x  by  virtue  of  "streaming"  in  seven-dimensional 
space,  (r,v,£) .  Since  the  time  rate  of  change  of  <J;  in  d'x  is  due  to  bubble 
streaming  into  d7x,  as  well  as  to  other  sources,  S,  we  have  3t|i/3t  =  F  +  S.  Using 
the  identity  V  •  (a^)  =  a-V  ty  +  ij»v  «a  and  the  fact  that  r  and  v  are  independent 
variables  (V  *v  =  o)  yields  equation  (7) . 

The  term  I  ty   in  the  bubble  transport  equation  is  generally  nonzero  due  to  the 
functional  dependence  of  acceleration,  a(r,v,£,t),  and  radius  change  rate, 
v(r,v,£,t),  on  the  bubble  velocity  and  radius,  respectively.   This  is  unusual  since 
in  transport  equations  normally  encountered,  the  term  corresponding  to  Z.ty   is 
generally  zero  (Chapman  &  Cowling  1964,  p. 46  and  p. 322) .  Mathematically,  E+  f  o 
causes  the  volume  element 

dV  =  [l+Etdt+0(dt2)]d7x 
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to  change  appreciably  in  the  time  interval  dt.  Physically,  the  term  -Z.ty   on 

the  right  hand  side  of  (7)  acts  as  a  pseudo  source,  or  sink,  depending  upon 

the  functional  form  of  a  and  v. 

For  example,  suppose  a  is  proportional  to  (-v)  .  Then  V  «a  <  o,  and  the 

bubbles  are  decelerated  into  a  smaller  volume,  d3v'  =  (1+V  «adt)d3v.  This 

increases  the  density  so  -(V  «a)  ip  acts  like  a  "source"  in  (7)  .  Similarly, 

suppose  v  is  proportional  to  v  I   so  that  8v/8£  is  proportional  to  v  .  For 

z  z 

v  >  o  the  bubble  radii  are  increasing  and  becoming  more  spread  out  over  a 

z 

larger  interval,  d£'  =  (l+(3v/3£)dt)d£.  This  decreases  the  density,  so  -(9v/9£)^ 
acts  like  a  "sink"  in  (7)  . 

The  second  term  on  the  right  hand  side  of  (7) ,  S(r,v,£,t) ,  represents 
distributed  bubble  sources  and  sinks.  In  addition  to  external  sources  and 
sinks,  this  term  can  include  the  effects  of  phenomena  such  as  bubble  scattering 
off  small  scale  turbulent  eddies,  bubble-bubble  interactions,  etc.  For 
instance,  when  the  scale  of  turbulent  eddies  is  the  order  of  bubble  radius,  they 
will  be  referred  to  as  "small  scale,"  and  bubbles  entrained  by  the  fluid  may 

scatter  off  of  them.  To  handle  this  phenomenon  one  could  define  a  turbulent 

scattering  cross  section,  Z  ,  such  that 

Zs(r,v1^v,£1^,t)v1^(r,v1,£1,t)d3v1d£1d7T 

3       * 
represents  the  probable  rate  at  which  bubbles  in  d  r  about  r  at  tune  t  are 

scattered  by  turbulent  eddies  from  d3v,  about  v,  and  d£,  about  i.   into  d3v 

about  v  and  dl   about  I,   where  v,  =  |v,|  is  the  speed.  Then 

S(r,v,£,t)  =   /  Z  (ZfV^fi^lf^v^irfVylyt)^^^^ 

(9) 
-v\|j(r,v,ji,t)   /  Z  (r^/W^tJd^d^ 
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represents  the  net  rate  at  which  bubbles  are  introduced  into  d7x.  Specification 
of  E  requires  detailed  knowledge  of  the  turbulent  field  as  well  as  a  model  for 
bubble-eddy  interaction. 

Substitution  of  equation  (9)  into  (7)  results  in  an  integro-differential 
bubble  transport  equation.  If  bubble-bubble  interactions  are  included,  the 
bubble  transport  equation  is  nonlinear  as  well  (Chapman  and  Cowling  1964,  p. 63) . 
Below  the  surface  in  the  upper  ocean  the  average  separation  between  bubbles  is 
about  1000  times  the  bubble  radii  (Medwin  1970)  ,  but  very  near  or  at  the  surface 
wave  action  can  cause  higher  bubble  densities.  In  addition,  the  energy  contained 
small  scale  turbulence  appears  to  increase  with  wave  action  but  to  decrease  with 
depth  (Shonting  1968)  .  To  avoid  at  this  time  obvious  analytical  complications, 
the  bubble  density  will  be  assumed  low  enough  so  that  bubble-bubble  interactions 
are  negligible.  It  will  also  be  assumed  that  the  energy  contained  in  small  scale 
turbulence  is  negligibly  small,  so  the  bubbles  are  partially  entrained  by  the 
fluid  without  abrupt  changes  in  velocity  or  radius.  These  assumptions  should  be 
valid  except  very  near  the  surface  in  heavy  seas. 

Suppose  the  source,  S(r,v,£,t),  is  independent  of  \\>   or,  at  most,  is 
proportional  to  \\>, 

S  =  Sl  +  V  *10) 

This  mathematical  simplification  yields  a  linear,  first  order,  partial  differ- 
ential bubble  transport  equation  whose  formal  solution  is  outlined  in  the  next 
section.  It  is  worthwhile  to  note  that  the  simplifying  assumption  (10)  is  not 
a  severe  limitation  for  many  problems  involving  gas  bubble  transport  in  the 
upper  ocean.  All  external  bubble  sources  such  as  surface  waves,  rain  and  snow 
nucleation,  organic  decay,  and  photosynthesis  can  be  included,  as  well  as  simple 
models  of  bubble  scattering,  absorption,  and  creation  that  have  the  form  (10)  . 
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III.   SOLUTION  BY  METHOD   OF  CHARACTERISTICS 

The  first  order  partial  differential  equation  (7)  is  equivalent  to  the 
following  set  of  eight  simultaneous,  first  order,  ordinary  differential 
equations  with  initial  conditions: 

dr/dt  =  v,  —  <  x,y,z  <  <*>,         r(t  )  =  rQ  =  (x  ,y  ,z  )      (lla-c) 

dv/dt  =  a,  -°°  <  v  ,v  ,v  <  oo       v(t  )  =  v  =  (v  ,v  ,v  )        (lld-f) 
x  y'  z  o    o    xo  yo  zo 

d£/dt  =  V,  O  <_  I    <   oo,  £,(t  )  =  l  (llg) 

dtp/dt  =  S-It*,  O  <  i>  <   »,  i|>(to)  =  4>Q  (12) 

The  equivalence  between  these  "characteristic  equations"  and  the  bubble  trans- 
port equation  is  most  lucidly  demonstrated  with  a  geometric  argument  given  by 
Garabedian  (1964,  p. 18) .  The  parametrically-represented  curve  C  in  nine- 
dimensional  space, 

S(t)  =    [t,x(t),y(t),z(t),v  (t),v  (t),v  (t),£(t),iMt)],  -oo  <  t  <  oo, 

obtained  by  integrating  the  characteristic  equations  (11)  and  (12) ,  is  called 
a  "characteristic  curve."  Components  of  the  initial  point,  R(t  ) ,  are  called 
"characteristics . "  The  family  of  characteristic  curves, 

{S(t):-»o  <  t  <  oo,  S(to)eP), 

arameterizedby  the  set  £  of  accessible  characteristics,  forms  a  hypersurface 
in  nine-dimensional  space,  ^(t;ro,vQ,Jlo,t  ),  that  satisfies  the  bubble  transport 
equation. 

The  family  of  characteristic  curves  is  essentially  parameterized  by  the  set 
(r  ,v  ,%  )   of  accessible  initial  conditions  since  t  and  ^0^0'v0'z0^  >  the 
initial  distribution,  are  generally  fixed  for  a  given  problem.  Since  a,  v, 
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and  S  are  generally  functions  of  r,v,£,  and  t,  the  solutions  to  (11)  and  (12) 
are  not  only  functions  of  the  parameter  t,  but  they  also  depend  upon  all  the 


initial  conditions  (r  ,v  ,1   ) : 

o  o  o 

r  =  r~(t;r"o,vo,Jlo,to)  e  r(t)  (13a) 

v  =  v(t;ro,vo,£Q/to)  =   v(t)  (13b) 

l  =   £(t;ro,vo,£o,to)  =  l(t)  (13c) 

and  i>  =   Mt;?^,^,^)  =  i/j  (t)  (13d) 

where  the  short  notation  on  the  right  is  used  for  brevity. 

The  characteristic  equations  (11)  are  simply  the  bubble  dynamics  equations. 
The  bubble  dynamics  model  used  to  construct  a  and  v  (see  $4)  ,  as  well  as  the 
characteristics  (r  ,v  ,1   ),  determine  a  "trajectory"  (r(t)  ,v(t)  ,£  (t) )  which  is 
the  projection  of  a  characteristic  curve  from  nine-dimensional  space  onto 
seven-dimensional  space.  Bubble  trajectories  in  (x,y,z)  space  are  likewise 
projections  of  characteristic  curves. 

A.  Lagrangian  and  Eulerian  Representations 

The  parametric  representation  (13) ,  obtained  by  integrating  the  character- 
istic equations  simultaneously,  describes  the  continuum  from  a  Lagrangian  frame 
of  reference.  In  this  description,  bubbles  are  labeled  by  group  according  to 
their  characteristics  (r  ,v  ,l   ) .  This  moving  frame  follows  a  specific  group, 

initially  in  d7x  =  d3r  d3v  d£  about  (r  ,v  ,%  )   at  t  ,  along  a  characteristic 
J  o  o   o  o        o  o  o     o     3 

curve  to  d7x  (t)  about  (r(t)  ,v(t)  ,£(t))  at  t.  At  any  time  t,  ip(t)d7T(t)  is  the 
number  of  bubbles,  belonging  to  the  group  labeled  (r  .v  ,£  ),  that  are  found 
ind7x(t)  about  (r  (t)  ,v(t)  ,1  (t) )  . 
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On  the  other  hand,  the  distribution  \p(r,v, l,t) ,  expressed  with  r,v,£,  and 
t  as  independent  variables,  more  conveniently  describes  the  continuum  from  an 
Eulerian  viewpoint.  The  Eulerian  frame  is  fixed,  and  attention  is  focused  on 
the  particular  volume  d7x  about  the  point  (r,v,l)  .  The  instantaneous  bubble 
density  iMr,v, £,t)dxT  is  considered  without  regard  to  the  characteristics  of 
bubbles  in  the  volume. 
Formal  Solution 


If  the  source,  S(r,v,£,t),  has  the  form  (10),  then  (7)  is  analogous  to  the 
Boltzmann  transport  equation  written  for  neutral  particles  (e.g.,  photons  or 
neutrons)  or  charged  particles  (e.g.,  electrons)  in  a  purely  absorbing  medium 
with  distributed  sources  (Case  and  Zweifel  1967,  p. 31  and  p. 247) .  For  example, 
neutrons  will  stream  along  their  characteristic  curves  in  (x,y,z)  space,  which 
are  straight  lines,  until  they  are  absorbed  or  until  they  escape  the  medium. 
Similarly,  the  bubbles  described  by  (7)  stream  along  their  characteristic 
curves  in  (r,v,  £)  space  until  they  disappear  (£  ->  o)  or  until  they  escape  the 
medium.  However,  the  characteristic  curves  defined  by  (11)  are  generally  not 
straight  lines  in  (r,v,il)  space  or  any  of  its  subspaces,  such  as  (x,y,z)  . 

With  S  given  by  (10) ,  the  formal  solution  to  (7)  is  obtained  by  integrating 
(12): 

+   /sl(?(t,),v(t•)£(t'),t•)exp[T1(t,t,)-T2(t,t,)]dt,    (14) 


where  V^V  ~     /  ^(f)  fv(f)  ,MO  ,t' )dt' 

to 

and  V^V   E  /  ^(t'Kvft'KJKfht'Jdt'. 
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All  integrals  in  this  expression,  and  all  similar  integrals  appearing  elsewhere 
in  this  paper,  are  understood  to  be  line  integrals  along  the  characteristic 
curve  that  passes  through  (r  ,v  ,%  )   and  (r,v,£) . 

To  obtain  the  number  of  bubbles  in  d7x  about  (r,v,Ji)  at  time  t,  (14)  is 
multiplied  by  this  volume  element.  However,  a  consistent  representation  for 
the  right  hand  side  of  the  equation  requires  use  of  the  Jacobian, 

Y(t,f)  e  3Cr{t)/v(t)/£(t))/3(r(tr),v(tl),Jl(t,)f), 

which  transforms  the  volume  element  from  one  point  to  another  on  the  same 
characteristic  curve: 

d7x  =  y^fJdV. 

Using  the  expression  derived  in  the  appendix, 

Y(t,t')  =  exp[T2(t,t')],  (15) 

the  solution  can  be  rewritten  in  the  convenient  form 

iMr,v,£,t)d7x  =  exp[T1(t,tQ)]iJ;(ro,vo,Jlo,to)d7To 

ft  (16) 

+  /  exptT^t'ns^r'^i^t'JdVdt', 

where  the  abbreviations  r  =  r(t;r  ,v  ,1   ,t  ) ,  r1  =  r(t';r  ,v  ,1   ,t  ) ,  etc.,  have 

o  o  o  o  o  o  o  o 

been  used.  From  this  expression  it  is  apparent  that  the  number  of  bubbles  in 
d7t  about  (r,v,£)  at  time  t  is  a  sum  along  the  characteristic  curve  that  passes 
through  (r,v,£)  at  t.  It  is  the  number  originally  in  d  t_  about  (r  rv _,*  )  at 
t  plus  (minus)  those  added  (subtracted)  by  sources  (sinks)  along  the  character- 
istic curve  between  (r  ,v  ,1  )   and  (r,v,£) .  If  E,  ^  o  in  (10)  then  the  result 
is  amplified  (attenuated)  by  the  factors  exp[T,]  because  of  the  creation 
(absorption)  rate  E,i|;  along  the  path.  For  t  =  t  +  dt,  (16)  expanded  in  a  Taylor 
series  to  0(dt)  yields  equation  (2) ,  as  it  must. 
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Equations  (7)  and  (16)  are  equivalent  statements  about  the  physical  model 
outlined  at  the  beginning  of  section  two.  They  formally  represent  relation- 
ships that  exist  between  the  bubble  distribution,  single  bubble  dynamics,  and 
the  bubble  sources.  As  such,  these  equations  can  be  used  for  calculations  to 
investigate  relationships  among  the  three  in  real  situations.  For  example, 
if  any  two  are  specified,  properties  of  the  third  can  be  inferred,  either 
directly  or  indirectly. 

Prior  to  using  (16)  for  calculations,  a  bubble  dynamics  model  (a  &  v)  must 
be  specified  so  that  equations  (11)  can  be  integrated  to  obtain  the  character- 
istic curves.  A  fair  amount  of  information  is  available  concerning  single 
bubble  dynamics,  and  in  the  next  section  some  of  the  important  points  are 
summarized . 

The  distributed  sources  and  fluid  properties,  as  well  as  the  appropriate 
bubble  dynamics  model,  are  dependent  on  the  particular  situation  being  analyzed. 
In  section  five  the  problem  of  gas  bubble  transport  in  the  upper  ocean  is  used 
as  the  basis  for  some  general  remarks  concerning  bubble  sources  and  fluid 
velocity  fields.  Then  equations  (7)  and  (11)  are  written  for  a  simple  ocean 
model,  and  some  calculations  are  performed  to  demonstrate  how  (16)  can  be  used 
to  evaluate  the  bubble  distribution  when  the  source  function  and  bubble 
dynamics  have  been  specific. 

It  is  significant  to  note  that  both  (14)  and  (16)  represent  the  distribution 
from  a  Lagrangian  view  point,  ^(t;r  ,v  ,£  ,t  ).  Transformation  to  the  more  con- 
venient  Eulerian  expression,  ij>(r,v,£,t) ,  requires  the  inverse  of  equations  (13), 
namely  r  (t  ;r,v,£,t),  v  (t  ;r,v,£,t),  and  l   (t  ;r,v,£,t).  In  general, 

O   O  O   O  O   o 

equations  (11)  will  have  to  be  integrated  numerically  to  establish  a  mapping 
between  (r  ,v  ,£  ,t  )  and  (r,v,£,t),  but  often  approximations  can  be  made  that 
will  yield  analytical  expressions. 
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If  S  is  a  more  complicated  function  of  ij>  than  (10)  ,  the  formal  solution  (16) 
represents  an  integral  equation  for  ty.     When  analytical  methods  aimed  at  solving 
(7)  or  (16)  fail,  the  Neumann  series  solution  to  the  integral  equation  can  be 
used  to  generate  numerical  results,  if  the  series  converges  fast  enough. 
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IV,   SINGLE  BUBBLE  DYNAMICS 

In  this  section  general  expressions  are  derived  for  a  and  v,  and  some  known 
models  of  gas  diffusion  and  drag  are  summarized  for  gas  bubbles  in  liquids. 
For  the  purpose  of  discussion  in  sections  four  and  five,  consider  an  inertial 
cartesian  coordinate  system  with  the  z-axis  vertical  upward  and  the,  origin  at 
the  surface  of  the  liquid.  Let  v  be  the  bubble  velocity  and  V(r,t)  be  the 
transporting  fluid's  velocity  relative  to  this  frame  of  reference. 
A.  Acceleration 

If  a  bubble  of  volume  ^containing  fluid  of  density  p,  were  completely 
entrained  (v  =  v)  by  a  fluid  of  density  p  ,  it  would  experience  the  same 
force,  p  adV/dt,  as  would  transporting  fluid  enclosed  in  the  same  volume.  For 
p  ^  p  ,  partial  entrainment  occurs,  and  the  bubble  maintains  a  relative 
velocity,  u  =  v-V,  with  respect  to  the  transporting  fluid.  Hence,  the  bubble 
experiences  a  drag  force,  F  ,  which  is  discussed  shortly  in  the  context  of 
gas  bubbles  in  liquids. 

As  the  bubble  moves  relative  to  the  transporting  fluid  it  "drags"  with  it 
an  amount  of  fluid  having  a  volume  equal  to  some  fraction,  3,  of  the  bubble 
volume.  By  Newton's  third  law,  this  adhering  mass  provides  a  supplemental 
reaction  term,  -  ( 3  a  p  du/dt) ,  that  acts  to  increase  the  bubble '  s  effective 
inertial  mass.  In  addition,  there  is  a  net  buoyant  force,  ag(p  -p)k,  where 
k  is  a  unit  vector  vertically  upward  and  g  is  the  acceleration  of  gravity. 
Combining  all  of  these  forces  in  Newton's  second  law  yields  the  following 
expression  for  acceleration: 


=  ?D/(apQ(3+p/po))  +  (l+3)(dV/dt)/(3+p/po) 


+  d-p/P0)gk/(3+p/P0)  (17) 
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in  which  dV/dt  =  8V/3t  +  (v-V  )^  (18) 

B.  Drag  on  Gas  Bubbles  in  a  Liquid 

Because  of  its  theoretical  and  practical  importance,  the  motion  of  gas 
bubbles  in  liquids  (e.g.,  air  bubbles  in  water)  has  been  actively  studied 
(Levich  1962,  4  80) .  The  regimes  of  bubble  motion  are  classified  according  to 
Reynold's  number,  R  =  up  SL/r\,   where  n  is  the  liquid's  viscosity.  Observations 
of  different-size  bubbles  rising  in  various  liquids  indicate  that  small 
(R  <  1)  and  medium-size  (1  <  R  <  700)  bubbles  maintain  a  spherical  shape, 
whereas  large  (R  >  700)  bubbles  deform  to  flattened  ellipsoids  and  very  large 
bubbles  (R  >  4500)  are  unstable  and  tend  to  break  up. 

For  small  bubbles  (I  <   lOOy  =  10_1+m  in  water)  the  drag  is  viscous  in  nature 
and  is  given  by 

-  Fv  =  -  <  n  l   u,  (19) 

with  k  =  4tt  for  a  perfectly  clean  bubble  having  a  mobile  two  fluid  interface 
at  its  surface.  However,  most  liquids  such  as  water  contain  "surface-active" 
materials  that  coat  the  bubble  and  destroy  the  mobility  of  this  interface.  In 
this  case  the  relative  velocity  of  the  transporting  fluid  goes  to  zero  at  the 
bubble's  surface,  and  the  bubble  behaves  like  a  solid  sphere  where  drag  is 
given  by  (19)  with  k  =   6tt  (Levich  1962,  ft  70  and  81)  . 

The  flow  past  medium-size  bubbles  (lOOy  <  I  <   2000y  in  water)  is  separated, 
with  the  separation  region  occupying  an  area  s^^  on  the  downstream  portion  of 
the  bubble's  surface.  Up  to  the  separation  point  resistance  that  acts  on  the 
bubble  is  viscous  in  nature  with  a  contribution  to  the  total  drag  given  by  (19) 
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with  k  =  12 tt.  Past  the  separation  point  and  into  the  bubble's  wake  the  flow 
is  characterized  by  turbulent  notion  with  a  contribution  to  the  total  drag 
given  by 

F"T  =  -O.SK^Sjiiu,  (20) 

where  K,.  is  a  drag  coefficient  with  0.65  >  K£  >  0.55  for  200  <  R  <  1000 
r  re 

(Levich  1962,  ^  80  and  82) . 

In  the  absence  of  surface  active  materials,  s,  is  very  small 

(s,  ^  *2/Re  for  R  >  >  1)  /  so  flow  past  the  bubble  is  essentially  unseparated 

and  the  drag  is  given  by  (19)  with  k  =  12-rr.  On  the  other  hand,  when  an  area 

s  on  the  bubble  is  covered  by  a  monolayer  of  surface-active  material,  the 

relative  fluid  velocity  in  this  region  is  zero,  and  flow  separation  occurs 

there.  In  this  case  viscous  drag  (19)  is  accompanied  by  the  form  drag  (20) 

with  s.,  =  s  ,  and  the  latter  dominates  the  former  when 
1    o 

so/(4tt£*)  >  28(n/pQ)2/(g*3) 

(Levich  1962,  $  82)  .  For  example,  when  I   =  500y,  only  2%  of  the  bubble's 
surface  need  by  covered  by  a  monolayer  of  this  material  before  the  drag  force 
becomes  quadratic  according  to  (20) . 

C.  Radius  Change  Rate 

For  gas  bubbles  moving  in  a  liquid  such  as  water,  radius  changes  are  caused 
primarily  by  changes  in  pressure,  due  to  changes  in  depth  or  surface  tension, 
or  by  gas  diffusion  across  the  bubble  surface.  Spherical  bubbles  (R  <  700) 
that  contain  n  moles  of  ideal  gas  at  temperature  T  and  pressure  P  are  described 
by 


4tt£3/3  =  nRT/P,  (21) 
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where  R  is  the  gas  constant.  For  isothermal  processes,  (21)  is  readily 
differentiated  to  yield 


u  =  l£  =  !  fdn/dt  _  dP/dt\ 

v  "  dt   3  I  n      P  / '  UZ) 


Gas  diffusion  is  included  in  the  dn/dt  term,  while  compression 
and  surface  tension  are  included  in  dP/dt.  Assuming  quasi-equilibrium,  gas 
pressure  is  the  sum  of  atmospheric  pressure,  P  ,  surface  tension  pressure,  2t,/l, 
and  fluid  weight,  -p  gz,  where  -z  >_  o  is  the  depth  and  c,   is  the  surface  tension: 

P  =  PQ  "Pogz  +  2?A,  (23) 


and 


dP/dt  =  -p  gv  -  (2cA2)d£/dt.  (24) 


D.  Gas  Diffusion 

As  results  of  the  next  section  will  indicate,  gas  diffusion  is  one  of  the 
more  important  factors  affecting  the  bubble  distribution.  It  is  a  complex 
phenomenon     depending  upon  many  influences  that  include  the  type  of  gas  in 
the  bubble,  the  gas  diffusivity,  D,  in  the  transporting  liquid,  the  gas  con- 
centration, C,  in  liquid  contacting  the  bubble's  surface,  the  gas  concentration, 
C  ,  in  liquid  far  away  from  the  bubble,  the  presence  of  surface-active  materials 
in  the  liquid  and  on  the  surface  of  the  bubble,  the  flow  field  around  the  bubble 
and  the  bubble  radius.  To  simplify  the  discussion  here  it  will  be  assumed  that 
gas  inside  the  bubble  is  composed  of  one  type  of  "average"  molecule  having  a 
single  diffusivity,  D. 

If  the  Peclet  number  is  large,  P  =  (u£/D)  >  >  1,  the  gas  concentration 

gradient  is  confined  to  a  thin  boundary  layer  on  the  surface  of  the  bubble,  and 

diffusion  between  C  and  C  occurs  across  a  small  distance  d  %  i/P  l/2.     Hence 

oo  e 

(Levich  1962,    tf  14  and  72) 
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dn/dt  =  -Gi2  (C-CJ  ,  (25) 

where  the  coefficient  G  is  generally  a  function  of  a,   u,  and  D.  For  specific 
applications  in  which  individual  components  of  gas  must  be  monitored,  equation 
(25)  will  have  to  be  generalized. 

Typically  (n/p)/D  ^  103,  so  there  is  a  wide  range  of  bubble  radii  for  which 
the  Peclet  number  is  large  (P  >  >  1)  but  the  Reynold's  number  is  small  (R  <  1) 
For  that  case  (small  bubbles) , 

G  =  8[U/6)  (DuA)]1/2  (26) 

in  the  absence  of  surface-active  material  (Levich  1962,  $$    72  and  91)  .  In  the 
presence  of  surface-active  material,  the  flow  field  around  a  small  bubble  is 
like  that  around  a  solid  sphere,  and  G  becomes  (Levich  1962,  s*  14) 

G  =  8(D2uA2)1/3   •  (27) 

A  representative  value  for  u  is  the  terminal  speed,  ll,  of  small  bubbles  rising 
in  a  quiescent  fluid.  Equating  the  drag  (19)  to  the  buoyant  force,  4ir£3gp  /3, 
which  assumes  p/p  <  <  1,  results  in  the  expression 

i^  =  (4tt/3<)  (gp0*2/n)  ,  (28) 

where  k  =  4-rr  in  the  absence  of  surface-active  materials,  and  k  =  6tt  in  their 
presence* 
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At  moderate  Reynold's  numbers  (medium-size  bubbles)  the  flow  is  separated, 
but  the  region  of  separation,  s, ,  is  generally  small.  If  gas  diffusion  across 
s,  can  be  neglected,  (25)  represents  the  gas  diffusion  with  G  equal  to  (26) 
multiplied  by  /3  (Levich  1962 ,  $  91) . 

If  liquid  in  contact  with  the  bubble  surface  is  assumed  to  always  be 
saturated  with  gas  at  the  pressure  P  (in  atmospheres)  prevailing  inside  the 
bubble,  then  C  =  KP,  where  K  is  the  absorption  coefficient  (kg*m~3'atm-1)  for 
the  liquid-gas  system.  Defining  the  partial  pressure  f  =  C^/K  (in  atm)  of  gas 
in  the  transporting  fluid,  (25)  can  be  rewritten 

dn/dt  =  -<2a2(P-f).  (29) 

Under  laboratory  conditions,  Wyman,  et  al  (1952) ,  measured  the  rate  of 
solution  and  the  change  in  composition  of  air  bubbles  in  stirred  sea  water,  as 
a  function  of  depth.  These  experimental  results  follow  very  closely 

dn/dt  =  -47r6£2(P-f)  (30) 

with  6=  4.4  x  10~5  moles  m~2 sec"1  atm-1.  Using  equations  (26)  and  (27)  with  pub- 
lished values  of  D  and  K  for  oxygen  and  nitrogen  in  water,  estimates  of  GK  are 
found  comparable  to  4tt6,  but  generally  larger.  The  presence  of  surface-active 
materials  that  retard  gas  diffusion  and  reduce  the  effective  absorption 
coefficient  across  the  bubble's  surface  may  have  contributed  to  the  discrepancy 
(Fox  &  Her zf eld  1954) . 
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V,   APPLICATION  TO  THE  UPPER  OCEAN 
The  purpose  of  this  section  is  to  briefly  discuss  application  of  bubble 
transport  theory  to  the  description  of  gas  bubble  populations  in  the  upper 
ocean.  Because  it  involves  a  variety  of  interrelated  phenomena  pertaining  to 
different  disciplines,  this  is  a  difficult  subject  to  analyze.  Also,  there 
is  only  limited  information  available  upon  which  to  base  mathematical  models 
of  bubble  sources,  fluid  velocity  fields,  gas  diffusion  rates,  etc.  However, 
it  is  the  author's  opinion  that  bubble  transport  theory  is  a  useful  tool  for 
interpreting  data  and  relating  phenomena  that  affect  the  distribution. 

To  begin,  it  is  appropriate  to  summarize  some  observations  concerning 
bubble  sources,  S,  and  fluid  velocity  fields,  V,  in  the  upper  ocean  and  to 
make  a  few  remarks  concerning  analysis.  Then,  to  demonstrate  an  application 
of  bubble  transport  theory,  some  sample  calculations  are  performed  and 
interpreted. 

A.  Bubble  Sources  and  Fluid  Velocity  Fields 

Physically  one  might  consider  bubble  sources  to  be  segregated  into  three 
categories  according  to  depth:  Those  on  the  ocean  floor,  those  concentrated 
at  the  surface,  and  those  distributed  throughout  the  medium,  principally  near 
the  surface.  Mathematically  these  classifications  are  convenient  for  analysis 
since  sources  in  the  first  two  categories  can  be  incorporated  into  (7)  and  (16) 
either  as  boundary  conditions  or  as  plane  sources. 

The  first  category  includes  sediment-initiated  bubbles  which  often  contain 
organic  gas.  For  example,  McCartney  and  Bary  (1965)  measured  relatively  large 
bubbles  with  radii  between  450y  and  800y  coming  off  the  bottom  of  Saanich  Inlet, 
British  Columbia,  and  speculated  that  these  bubbles  contained  methane. 
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In  the  second  category  are  bubbles  generated  at  or  very  near  the  surface 
by  wind  and  wave  action,  precipitation,  and  captured  aerosols.  The  subsurface 
measurements  of  Blanchard  and  Woodcock  (1957) ,  Glotov,  et  al  (1962) ,  and 
Medwin  (1970)  imply  that  most  of  the  bubbles  associated  with  breaking  waves 
have  radii  less  than  200y.  On  the  other  hand,  experiments  conducted  by  Monahan 
and  Zietlow  (1969) ,  in  which  a  small  volume  of  seawater  was  poured  into  a  tank 
of  the  same,  indicate  that  the  cloud  of  bubbles  formed  by  the  descending  plume 
of  mixed  air  and  water  contains  an  appreciable  number  of  200y  to  1000  y  bubbles 
as  well.  In  addition  to  wave-generated  bubbles,  rain  and  snow  nucleate  bubbles 
at  the  surface  which  have  radii  predominantly  less  than  50y  (Blanchard  and 
Woodcock  1957) .  To  explain  the  excessive  sound  attenuation  and  scatter 
observed  in  near-surface  water  at  high  frequencies  (f  >  60kHz) ,  Medwin  (1970) 
speculated  that  continental  aerosols  may  be  the  primary  source  of  small  bubbles 
U  <  50/u)  near  the  surface. 

In  the  third  category  are  bubbles  initiated  within  the  volume  of  the  upper 
ocean.  Internal  wave  action,  radiation,  and  biological  activity,  such  as  photo- 
synthesis and  marine  life  gas  emissions,  have  been  postulated  to  be  responsible 
for  many  of  the  smallest  bubbles  (£  <   60 y)  .  Also,  some  of  the  volume  sources 
may  be  associated  with  pressure  fluctuations  caused  by  wind-driven  surface 
waves.  Pressure  changes  and  turbulence  may  nucleate  small  bubbles  at  distri- 
buted cavitation  sites,  while  the  vertical  fluid  velocity  and  turbulence 
created  by  breaking  waves  introduces  larger  bubbles  below  the  surface. 

When  a  wave  breaks  and  air  is  rapidly  mixed  with  water,  the  bubbles  formed 
are  entrained  by  a  strong  vertical  velocity  field;  these  bubbles  have  been 
observed  as  far  as  three  wave  heights  below  the  surface  (Kanwisher  1963)  .  As 
the  strong  downward  currents  decay,  the  bubbles  tend  to  rise  under  buoyancy 
with  motion  modified  by  the  remaining  fluid  velocity  field  (ref .  equation  (17))  . 
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The  largest  bubbles  rise  most  quickly  and  reach  the  surface  before  going  into 
solution  while  the  smallest  bubbles  succumb  to  gas  diffusion  (LeBlond  1969) . 
As  a  result,  the  most  prevalent  bubble  within  a  few  meters  of  the  surface  should 
have  a  radius  between  the  two  extremes,  an  assertion  supported  by  the  data  of 
Medwin  (1970)  and  Glotov,  et  al  (1962)  who  measure  distribution  maxima  at  90y 
and  60y,  respectively. 

Shonting  (1968)  has  measured  the  autospectra  of  particle  motions  in  the 
upper  ocean  and  shown  that  the  dominant  peaks  are  associated  with  the  fre- 
quencies of  the  surface  waves.  The  spectral  energy  decays  exponentially  with 
depth  in  a  way  that  attenuates  high  frequencies  the  most  rapidly. 

Example 


Consider  a  one-dimensional,  steady  state  ocean  in  which  all  velocities  are 
vertical  and  in  which  all  functions  depend  spatially  on  depth  only.  In  this 
case  the  bubble  transport  equation  reduces  to 

where ^  (z,v,£)  is  the  distribution  function,  S(z,v,£)  is  the  distributed  source, 
v  =  v  is  the  vertical  velocity,  a  =  a(z,v, I)   is  the  vertical  acceleration, 
v(z,v,ft)  is  the  radius  change  rate,  and  I.   =  3a/3v  +  3v/3£. 

The  characteristic  equations  for  (31)  can  be  written  with  z  <^  o  chosen  as 
the  independent  parameter.  This  choice  is  equivalent  to  dividing  equations 
(11)  and  (12)  by  (lie)  : 

dv/dz  =  a(z,v,£)/v,       -00  <  v  <  °°     &    V(ZQ)  =  v0;  ^32a^ 

dJl/dz  =  v(z,v,£)/v,       o  <_  l  <  «  &  Mz  )   =  *Q;  (32b> 

d<f>/dz  +  (£t/v)<j>  =  S(z,y,£)/V/  o  <_  <j>  <  °°  &  <fr  (zq)  =  4>0(v0'£0)    (33) 
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Simultaneous  integration  of  equations  (32)  yields  the  family  of  characteristic 
curves 

{z,  v(z;zo/vq/£o)/  £(z;zo/vq,<Io);  z  <^  o}        (34) 

and  bubble  trajectories  in  (x,y,z)  space  are  vertical  straight  lines. 

A  formal  solution  analogous  to  (14)  is  readily  obtained  by  integrating  (33) 
along  characteristic  curves  (34) .  The  form  analogous  to  (16) , 

v<}>(z,v,£)dvd£  =  v04>(z0>v0'£0)dv0d£0  +  /  S(z,,vl/£,)dzldvld5-'/   (35) 


is  derived  Using  the  Jacobian, 


Zo 


y(z,z')    =  l^fjpy  =  £-  exp   [    /    Zt(z\v",£")dz"/v"].      (36) 

Jz' 

Expression   (36)   follows  from  the  differential  equation 

dy/dz  =   (Zt  -    (a/v2))y  (37) 

obtained  when  characteristic  equations  (32)  are  combined  in  a  manner . similar  to 

the  general  scheme  detailed  in  the  appendix.  Of  course,  all  integrals 

appearing  in  this  section  are  taken  along  characteristic  curves  (34) ,  where 

abbreviated  notation  v'  =  v(z';z  ,v ,1  )      I"   =  £(z";z  ,v  ,l  ) ,   etc.  has  been 

used  in  (35)  and  (36)  . 

Physically  equation  (35)  states  that  the  flux  of  bubbles ,   v<j>,  at  the  point 

(z,v,£) , equals  the  flux  of  bubbles,  v  <J>  ,   at  the  point  (z  ,v  ,1   )  ,   plus  (minus) 

the  bubbles  added  (removed)  by  distributed  sources  (sinks)  along  characteristic 

curve  (34) .  Further  physical  insight  can  be  gained  by  reflecting  on  a  situation 

where  there  are  no  distributed  sources,  S  =  o,  so  that  <i(z;z  »v  I  )dvdl  = 

o  o  o 

(v  /v)d)(z  ,v  ,£  )dv  ds,  .  If  there  is  no  fluid  velocity,  bubbles  rise  under 
o     o  o  o   o  o 

buoyancy  at  terminal  velocity  and  may  shrink  and  slow  down  or  expand  and  speed 
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up  (LeBland  1969).  Those  belonging  to  the  Lagrangian  group  (z  ,v  ,l   )  whose 
velocity  and  radii  increase  with  z  have  a  density  that  decreases  as  z  increases 
because  the  bubbles  at  z  +  dz  are  running  away  from  those  at  z.  On  the  other 
hand,  those  whose  velocity  and  radii  decrease  as  z  increases  have  a  density 
increasing  with  z  since  the  bubbles  at  z  are  piling  up  on  those  at  z  +  dz. 

1.  Characteristic  Curves 

A  first  step  toward  using  (35)  to  calculate  the  bubble  distribution  is  the 
specification  of  a(z,v,£)  and  v(z,v,£)  and  the  computation  of  characteristic 
curves  (34) .  A  useful  expression  for  the  radius  change  rate  is  obtained  by 
substituting  (21) ,  (23) ,  (24) ,  and  (29)  into  (22)  and  solving  for  d£/dt.  Then 
(32b)  becomes 

d£/dz  =  v/v  =  {U/3)-(GK/4tt)  (KT)  [D(l-f)-z  +  (T/i)  ]/v}/[D-z  +  (2/3)  (r/£)]   (38) 

where  (-z)  is  the  depth,  f  is  the  partial  pressure  (in  atmospheres)  of  gas  in 
the  water, 

D  e  p  /p  g  =  10  meters 
a  Ko^ 

is  the  depth  of  sea  water  equivalent  to  one  atmosphere,  and 

r  e  2c/p  g  =  1.47xl0"5  meters2 

is  a  surface  tension  coefficient  in  water.  Note  that  the  T/l   terms  are 
negligible  unless  the  bubbles  are  very  small  (l   <  30y) . 
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Acceleration  is  given  by  (17)  with  the  drag  force  given  by  (19)  or  (20) . 
If  the  fluid  velocity  is  negligible  (V  =  o) ,  the  bubbles  rise  very  nearly  at 
their  terminal  velocity,  v  .  This  is  seen  in  Figure  1  which  was  obtained  by 
simultaneously  integrating  equations  (32)  numerically*  for  bubbles  of  various 
initial  radii,  I  ,   starting  from  rest  at  z  =  -20  meters.  Small  bubbles 
(£  <  lOOy) ,  and  medium-size  bubbles  (lOOy  <  I  <   2000y)  with  viscous  drag  domi- 
nant, rise  at  terminal  velocity  y _1  given  approximately  by  equating  (19)  to 
the  buoyant  force: 

v^  =  (g/(<'a)H2  (39) 

where  k1  =  1  and  k*  =  2,  respectively,  and  where 


a  =    (9/2) (n/p  )  =  4.37x10  6  meters2/sec 
in  water.  When  medium-size  bubbles  have  a  sufficient  fraction  of  their  surface, 

e  =   so/(4tt£2), 

covered  by  a  monolayer  of  surface  active  material,  turbulent  drag  dominates 
viscous  drag  because  flow  separation  occurs  at  the  boundary  of  s  .  They  rise 
at  terminal  velocity  vT?  given  approximately  by  equating  (20)  to  the  buoyant 
force : 

vT2  =  (g/O1/2*1/2  (40) 


where  5  5  (3/2) Kf  e  =  0.9e. 


*A11  numerical  integrations  alluded  to  in  this  section  were  carried  out  using 
gaussian  quadrature  or  optimal  fourth  rank  Runge-Kutta  (Ceschino  and  Kuntzmann 
1966,  p.  67). 
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Datta,  et  al,  (1950)  has  averaged  the  results  of  several  different  observers 
to  produce  an  empirical  curve  for  the  terminal  velocity  of  air  bubbles  in  water 
as  a  function  of  their  radii.  Values  produced  by  (40)  agree     with  this 
empirical  data  for  300y  <  i  <   lOOOy  if  the  surface  fraction,  e,  is  reduced  mono- 
tonically  as  the  radius,  I,   increases.  This  leads  one  to  speculate  that  as  the 
radius  and  terminal  velocity  increase,  rapid  liquid  motion  at  the  surface  either 
washes  away  surface  active  molecules  or  carries  them  to  the  rear  where  they  form 
a  compact,  saturated  monolayer  (Levich  1962,  p. 447) . 

In  general  the  terminal  velocity  of  a  rising  bubble  is  not  given  simply  by 
(39)  or  (40)  during  its  entire  lifetime.  Many  bubbles  shrink  or  expand  as  they 
rise  (LeBlond  1969)  and  thereby  transition  from  the  regime  of  predominantly 
turbulent  drag  to  the  regime  of  predominantly  viscous  drag,  or  vice-versa.  The 
treatment  of  separated  flow  in  Levich  (1962,  $  82)  suggests  that  the  sum  of  (19) 
and  (20)  should  be  a  useful  model  for  the  drag  on  medium-size  bubbles  in  the 
transition  region.  Using  this  model  in  (17) ,  (32a)  can  be  written 

dv/dz  =  -  [  (5/6)  A]  [  (v-v_)  (v+v+)  /v] ,  (41) 


where 


v+  =   v  [±1  +  (l+4Q2)1/2]/(2Q)  (42) 


with 


and 


Q  e  vT1/vT2  =  UA)3/2 


i  =  [(<-a)2/(g?)]1/3 

Numerical  integration  of  (38)  and  (41)  demonstrates  that  the  bubbles  released 
from  rest  rise  with  terminal  velocity  vT  =  v_  (figure  1) .  This  approach  appears 
to  have  some  merit  since  (42)  for  v_  yields  values  close  to  the  empirical  curve 
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by  Datta,  et  al  (1950)  and  since  the  model  shews  the  proper  limiting  behavior: 
For  £  <  I,   viscous  drag  dominates  turbulent  drag  in  (41)  and 

v_  =  vT]_  -  vT1{Q2  +  0(QI+)} 

for  I  >  I,   turbulent  drag  dominates  in  (41)  and 

v_  =  vT2  +  vT2{0.5QTl  +  0(Q"2) }. 

Computation  of  characteristic  curves  (34)  in  a  quiescent  fluid  (V  =  o)  is 
simplified  because  the  bubble  velocity  is  always  very  nearly  the  terminal 
velocity,  v=v  (z;z  ,1   ),  and  is  essentially  independent  of  v  .  Either  (39), 
(40) ,  (42),  or  empirical  values  for  v  can  be  used  in  (38)  to  obtain  £(z;z  ,£  ) 
If  (39)  or  (40)  is  used  and  surface  tension  is  neglected,  closed  analytical 
expressions  result:  For  v  =  v_, ,  r  =  o,  and  f  constant, 


,D-z 

(43a) 


-K'V-  {tr)*o3  +  j  W^[(z2-zo2)-2Da-f){z-zo^1/3' 


with  an  inverse  transformation 


£o  =  £Tl(zo;Z'£);  (43b) 


for  v  =  v  2,   r  =  o  and  f  constant, 


P~z0\  \/l  2./3 

-3D(2-f)  [l-lg-^j        ]]}  (44a) 
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with  an  inverse  transformation 

£o  =  AT2(zo;z'£)-  (44b> 

Otherwise,  (38)  must  be  integrated  numerically. 

To  evaluate  accuracy  of  approximate  analytical  results  (43)  and  (44) ,  l   n 
and  I   ~  are  plotted  vs  depth  in  figures  2  ,3,  and  4  and  compared  to  I    (z;z  ,1    )  , 
obtained  by  numerically  integrating  (38)  with  v  =  v_  given  by  (42)  .  For  the 
sake  of  example,  f  =  1  and  KTGK/4u  =  106  meters/sec  (Wyman,  et  al  1952)  are 
used  in  all  three  models,  and  £  ..  (z),  £T2(Z)  and  £T(Z)  are  represented  for 
various  initial  radii,  £  ,  at  z  =  -20  meters.  Note  that  I   ,  is  a  good 
approximation  to  a     for  30y  <  l  <   lOOy  (viscous  regime)  and  that  l  ~   closely 
approximates  £„  for  I  >   400y  (turbulent  regime) .  Neither  (43)  nor  (44)  is 
valid  for  I  <   30y  because  they  neglect  surface  tension,  and  neither  is  accurate 
in  the  transition  region  lOOy  <  %  <   400y. 

An  example  of  the  phenomenon  analyzed  by  LeBlond  is  evident  in  Figures 

2-4.    As  a  bubble  rises  under  buoyancy  it  tends  to  expand  as  the  hydrostatic 

pressure  decreases,  but  it  also  tends  to  contract  as  gas  diffuses  across  the 

bubble's  surface.  If  a  bubble's  radius  is  less  than  some  critical  radius,  L  , 

it  rises  too  slowly,  shrinks  monotonically,  and  disappears  before  reaching  the 

surface  because  gas  diffusion  dominates  hydrostatic  expansion.  If  its  radius 

is  greater  than  another  critical  radius,  £  ,  it  will  expand  monotonically  as 

a 

it  rises  because  hydrostatic  expansion  dominates  gas  diffusion.  For  an  initial 
radius  between  l     and  SL,   the  bubble's  behavior  is  not  necessarily  monotonic; 
it  may  or  may  not  reach  the  surface. 

Bubbles  with  radii  less  than  L    will  rise  a  finite  distance  before 
disappearing,  and  computed  values  of  I     and  JL  depend  upon  the  gas  diffusion 
TDdel  used  as  well  as  parameters  like  z  .  Examination  of  ^T  in  Figures  2-4 
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indicates  that  I     ^   400y  and  L    ^   250y,  while  a  lOOy  bubble  rises  about  1.1 
meters.  Unfortunately  there  is  usually  insufficient  information  available 
for  an  accurate  specification  of  the  gas  diffusion  phenomenon.  For  example, 
in  (29)  the  gas  diffusion  rate  depends  not  only  upon  the  partial  pressure  f (z) , 
but  also  upon  the  parameters  G  and  K  which  are  determined  by  the  type  of  gas 
and  the  presence  of  surface  active  materials.  This  data  is  often 
unavailable  for  the  ocean.  Since  characteristic  curves  and  bubble  distributions 
(ref .  Figures  5 -12) are  sensitive  to  the  gas  diffusion  rate,  specification  of 
the  unknown  parameters  is  important,  and  it  may  require  a  sequence  of  educated 
guesses  refined  from  comparisons  of  computed  and  observed  distributions. 
If  the  fluid  velocity  V  is  not  negligible  and  the  problem  is  multi- 
dimensional, the  computation  of  characteristic  curves  is  more  difficult.  In 
general,  equations  for  x,y,v  ,v  ,  and  v  must  be  integrated  simultaneously 
with  (38) .  However,  if  dV/dt  given  by  (18)  is  small  compared  to  the  other  terms 
in  (17) ,  a  useful  approximation  (that  eliminates  the  three  characteristic 
equations  for  velocity)  might  be 

v  ~  V  +  v^jc, 

where  v  is  the  terminal  velocity  of  the  rising  bubble. 

2.  Bubble  Distribution 

Since  bubbles  rise  in  a  quiescent  fluid  at  very  nearly  their  terminal 

velocity, 

4>(2,v,Jl)  =  <5(v-vT)  §(z,£)  (45) 
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should  be  a  good  approximation  where  vT  is  given  either  empirically  or  by  (39) , 
(40)  or  (42) .  In  this  equation 

§(z,£)  e  /   cf>(z,v,£)dv  (46) 

J     00 

is  the  "radius  density,"  the  number  of  bubbles  per  unit  volume  per  unit  radius. 
Equation  (35)  can  be  integrated  over  v  to  yield 

$(z,£)  =  (v^/v^  £(zo,£o)(3S,o/8£)  +  (l/vT)  J    s(z',£')  O^'/aiOdz'   (47) 

o 
with  « 

s(z,,£I)  ~j     S(z'/v,/£')dv,/  (48) 

■J  —00 

where  d£Q  =  (8£ /3£)d£  and  d£*  =  (3£'/3£)d£  have  been  used.  The  bubble  density, 

f" 

$(z)  =   /  §(z,£)d£,  (49) 

Jo 

which  is  the  total  number  of  bubbles  per  unit  volume  at  depth  z,   regardless  of 
radius,  is  given  by  (47)  integrated  over  all  radii, 

i(z)    =   I   (Vrib/vT)|(zo,£o)d£o  +   /    /  (l/vT)s(z',£')0£'/3ydz'd£o   . 

Jo  J  o  J  z 

°  (50) 

It  must  be  remembered  that  vm  and  V™_ ,  as  well  as  £  and  £ ' ,  are  functions  of  £ 

T     To  '  o 

and  that  integrations  in  z  are  taken  along  the  characteristic  curves 

{zI,vT(z,;zo,£o),£T(z';zo,£o):zo  <_  z'<_  z},         (51) 
where  £  is  obtained  by  integrating  (38)  with  v  =  v  . 


35 


For  most  applications  the  radius  density  provides  adequate  information. 
To  calculate  this  function  for  a  given  z  and  I,   one  needs  9£'/8£,  V  ,   and 
s{z\V)   at  points  on  the  characteristic  curve  that  passes  through  (z,l)   so 
the  integral  in  (47)  can  be  evaluated.  Using  notation  from  the  appendix, 

9£'/3£  =  1/y77(z,z'),  (52) 

where  y77  satisfies  the  differential  equation  obtained  by  setting  v  =  v_  in 
(38)  and  taking  the  partial  derivative  with  respect  to  V : 

dy77/dz  =  F?7  y77  (53) 

with  the  boundary  condition 


Y7?(z/z)  =  1 


where 


F77<z;zo,£o)  =  {(2r/£2)[£-(36RT/vT)(D(l-f)-z+rA)][3(D-z)  +  2T/l]~l  + 


[1+  (36RT/vT)  [T  /  Sl2+ (D  (1- f ) -z+T  / 1)  (Zv^/d  I) /vT]]}  [3  (D-z)   +  2T/i]^    (54) 

in  which  (GK/4tt)  =  constant  e  6  (Wyman,  et  al,  1952)  has  been  used.  When  the 
approximation  vT  =  v_  (42)  is  used, 

dv^/dSl  =   (3vTl/£)  [1  +  4Q2]-l/2  _  v_/£.  (55) 

The  solution  to  (53)  is 

Y77(z,z')  =  exp[7(z,z')]  (56) 


with 


J(z,z')    /  F77(z";zo,£o)dz", 
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and  it  is  dbvious  that 

Y77(z/zI)  =  y77(z,zo)/y77(z',zo).  (57) 

Using  (56)  and  (57) ,  (47)  can  be  written 

§(z,Z)   =  [exp[-7(z,zo)]A'T]{vTo$(zo,Jlo)  + 

Tz  (58) 

+  /   s(z,,£,)exp[y(z,,zo)]dz,> 

o 

Once  the  source  s(z,£)  has  been  specified,  (56)  can  be  used  to  evaluate  the 
radius  density  at  points  along  any  characteristic  curve  (51)  obtained  by 
integrating  (38) .  In  the  following  discussion  of  sample  calculations, 
$(z,£)  will  refer  to  the  radius  density  when  the  approximation  v  =  v  (42) 
is  used. 

If  either  approximation  (43)  or  (44)  is  appropriate,  then  either 

3&'/3«,  -  [(D-z)/(D-z')]  U2A'2]  (43c) 

or 

3A'/3Jl  =  [(I>-z)/(D-z,)]l/2UA,]1/2,  (44c) 

is  used  in  (47)  to  obtain  the  approximate  radius  density  $..(z,&)  or  §2(z,£), 
respectively,  at  points  along  the  characteristic  curve  (51)  with  &  =  £  , ,  or 
£T  =  I   ~,  respectively.  If  gas  diffusion  can  also  be  neglected  (GK  =  o)  then 
both  approximations  reduce  to 

I   =  £T3(z;zo,Jlo)  =  [(D-z^/fD-z)]1/3^  (59) 

which  can  be  used  in  (47)  and  (51)  to  obtain  the  approximate  radius  density 
$3(z,Ji). 


37 


3.  Numerical  Results 

For  the  sake  of  example,  consider  the  situation  investigated  by  McCartney 
and  Bary  (1965)  who  acoustically  measured  rather  large  gas  bubbles  ascending 
from  the  bottom  of  Saanich  Inlet,  B.  C,  and  inferred  the  bottom  source  repre- 
sented in  Figure  5  from  the  ascent  velocities.  If  we  model  that  ocean  as  1-D 
and  quiescent,  then  equation  (47)  can  be  used  to  infer  the  radius  density 
throughout  the  volume.  Since  volume  sources  produce  mainly  smaller  bubbles 
that  their  measurements  do  not  include,  only  the  plane  source  $  (z  ,1  ) 
(Figure  5 )  at  z  =  -197  meters  is  considered  in  this  sample  calculation,  and 
(47)  reduces  to 

$(z,£)  =  (vTb/vT)i(zo,£o)OJlo/35,).  (60) 

Using  f  =  1,  approximations  §,$,,$-/  and  <£_  to  the  radius  were  computed 
at  various  depths  between  z  =  -197m  and  z  =  0.   $ (z,£)  (solid  line)  and 
§L(z,£)  (dashed  line)  are  plotted  in  Figures  6-10   at  the  depths  z  =  -172m,~147r 
-122m,  -97m,  -47m,  and  Om.  Although  actual  calculations  do  not  yield 
precisely  rectangular  distributions  for  z  /  z  ,  $  and  $-  and  drawn 
as  perfect  histograms  because  the  error  introduced  is  negligible. 
While  £o  neglects  gas  diffusion  and  surface  tension,  J  includes  these 
effects.  For  the  gas  diffusion  model  used,  all  bubble  radii  for  which 
|(z  ,SL   )  is  nonzero  lie  below  the  critical  radius  L,   so  these  bubbles  will 
shrink  monotonically  as  they  rise  and  will  disappear  before  reaching  the  surface. 
The  distribution  $  is  seen  to  shift  toward  smaller  radii  as  depth  decreases, 
until  even  the  largest  bubbles  have  disappeared  by  the  time  z  ^  -90m,  so  only 
§^  is  shown  in  Figure  10.  since  JL  =  o  in  the  model  used  to  calculate  §.,  these 
bubbles  necessarily  expand  as  they  rise,  and  the  distribution  shifts  monotonically 
to  larger  radii. 
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The  bubble  density  is  simply  the  area  under  the  histogram.  As  z  increases 
from  -197m  the  density  of  $(z,Jl)  initially  increases  since  flux  is  conserved 
and  the  bubbles  are  slowing  down  and  piling  up.  However,  nearer  the  surface, 
most  smaller  bubbles  have  succumbed  to  gas  diffusion  and  the  density  starts  to 
decrease  toward  zero.  For  $_(z,£)  the  density  is  seen  to  decrease  as  z 
increases  since  the  bubbles  are  accelerating  and  becoming  spread  out. 

The  distributions  $2  and  *_  are  not  shown  because  neither  follows  $  closely. 

$2  is  the  best  approximation,  but  these  bubbles  do  not  quite  shrink  fast  enough. 

Some  of  the  bubbles  in  $,  contract  monotonically  while  others  contract  and  then 

expand.  For  the  model  used  most  of  the  bubbles  originate  at  the  bottom  with 

radii  in  the  transition  region,  I     <  l     <  L,   where  neither  approximations  (43) 

a   o   d 

nor  (44)  are  valid. 

Interestingly  enough,  the  experimental  results  of  McCartney  and  Bary  (1965) 
indicate  the  bubbles  from  the  bottom  of  Saanich  Inlet  expand  as  they  rise 
according  to  the  distribution  $_,  which  ignores  gas  diffusion.  Thus,  it  appears 
that  diffusion  is  somehow  inhibited  across  the  surfaces  of  these  bubbles  which 
are  composed  of  "gas  released  from  the  highly  organic,  anaerobic  sediment  of  the 
bottom, "  possibly  methane .  Either  the  water  was  saturated  with  respect  to  the 
gas (es)  in  the  bubbles,  or  these  bubbles  were  coated  with  an  active  material  that 
severely  impeded  gas  diffusion.  While  it  is  significant  to  note  that  $  is  cal- 
culated using  GK/4tt  =  6  (eqn.  (30)) ,  which  is  based  upon  experiments  with 
relatively  clean  air  bubbles  (Wyman,  et  al  1952) ,  and  that  f  =  1  at  all  depths 
is  not  physically  realistic,  certain  qualitative  conclusions,  based  on  a  com- 
parison of  $  and  $~  to  the  observed  distribution,  are  still  valid.  For  example, 
gas  diffusion  appears  to  be  negligible  in  this  situation;  and  gas  diffusion 
WDuld  profoundly  affect  the  bubble  distribution.  The  latter  observation  emphasizes 
the  fact  that  results  are  strongly  dependent  on  the  bubble  dynamics  and  ocean 

models  used. 
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As  the  basis  for  a  second  sample  calculation,  consider  the  supposition 
that  the  observed  bubble  distribution  indicates  the  functional  dependence  of 
the  volume  source.  To  test  this  idea,  consider  a  hypothetical  volume  source 
suggested  by  the  measurements  of  Shulkin  (1968)  and  Medwin  (1970) : 

s(z,£)  =  ez/hL(£)  (61) 

where 

(£/£')~3,5for  £  >  if 


L(£)  =  i 


(62) 
0,for  0  <_  £  <   £' 


with  h  =  10  meters  and  £'  =  20y.  Physically,  a  source  decaying  exponentially 
with  depth  might  arise  from  photosynthesis,  light  sensitive  bacteria,  wave 
action,  etc.  Medwin  (1970)  has  observed  an  exponential  behavior  for  small 
bubbles  (£  <  60y)  ,  and  both  Medwin  (1970)  and  Shulkin  have  observed  that  the 

near  surface  bubble  density  appears  to  fall  off  like  same  power  of  the  bubble 

-3  5 
radius.  The  radius  spectrum  L(£)^£  *  is  chosen  as  an  example. 

Assuming  that  the  source  is  zero  below  z  =  -30  meters,  and  that  f  =  1, 
the  radius  density  was  computed  using  v  =  v_  in  (38)  and  (58)  .  Figures  11  &  12 
sents  log-log  plots  of  $(z,£)  (solid  line)  and  s(z,£)  (dashed  line)  vs  radius 
for  four  depths,  z  =  -29m,  -20m,  -10m,  and  0m.  A  comparison  of  §  and  s  indi- 
cates that  the  distribution  follows  approximately  the  functional  form  of  the 
source  for  £  >  lOOy: 

I (z,£)  *  ez/llV3*5  (63) 

h1  =  6.5m.  The  distribution  at  z  =  -29  meters  does  not  fit  this  pattern 
because  it  is  too  near  the  region  assumed  to  be  devoid  of  bubbles. 

For  £  <  100y  the  functional  form  of  $  is  more  complicated  and  does  not 
necessarily  follow  (63) .  Except  at  the  surface,  the  curves  are  all  concave 
down  which  indicates  a  deficiency  of  small  bubbles.  Gas  diffusion  and  surface 
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tension  cause  small  bubbles  to  disappear  much  more  rapidly  than  large  bubbles. 
While  the  small  bubbles  observed  at  a  given  depth  originate  primarily  from 
local  sources,  many  of  the  larger  bubbles  observed  originate  at  sources  distri- 
buted well  below  that  depth.  However,  for  the  model  used,  gas  diffusion 
decreases  rapidly  with  depth,  so  there  is  less  discrimination  against  the  small 
bubbles  in  the  near  surface  region.  From  Figures  11  and  12, 

JtO#A)~  I   "3*75  (64) 

These  results  indicate  that  if  gas  diffusion  is  significant,  the  small 
bubble  sources  will  probably  fall  off  to  a  higher  power  in  I   than  the  observed 
distribution,  but  the  larger  bubble  sources  may  have  roughly  the  same 
£ -dependence  as  the  observed  distribution.  Furthermore,  the  depth  dependence 
of  the  bubble  sources  is  likely  to  be  functionally  similar  to  that  of  the 
bubble  population.  Of  course  these  qualitative  relationships  between  bubble 
population  and  bubble  source  are  not  likely  to  be  applicable  unless  the  model 
is  valid.  For  example,  any  of  the  following  can  have  a  profound  influence  on 
the  bubble  distribution  and  its  relationship  to  the  source:  A  significant 
fluid  velocity,  V  ^  o,  a  partial  pressure,  f (z) ,  that  varies  with  depth,  or 
surface  active  materials  that  impede  gas  diffusion  (e.g.,  lower  K  in 
equation  (24) ) . 

The  distributions  $,  and  <$2  were  also  computed  for  this  example  but  are 
not  displayed.  As  expected,  §,  approximates  $  fairly  well  for  30  y  <  a  <   lOOy, 
and  $2  is  a  good  approximation  to  $  for  a  >   350y. 
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VI.   CONCLUSIONS 
The  main  purpose  of  this  paper  has  been  to  develop  the  general  framework  of 
bubble  transport  theory  and  then  to  demonstrate  its  application.  The  bubble 
transport  equation  and  its  solution  relate  the  ensemble  average  behavior  of 
a  bubble  population  to  dynamics  of  a  single  bubble  and  to  properties  of  the 
transporting  fluid.  Among  the  most  important  aspects  of  gas  bubble  dynamics 
in  liquids  are  bubble  gas  diffusion  and  drag,  which  depend  upon  many  different 
physical  parameters.  For  the  upper  ocean,  information  concerning  these  para- 
meters, as  well  as  information  about  bubble  sources  and  fluid  velocity  fields, 
is  usually  incomplete.  In  this  case,  bubble  transport  theory  cannot  be 
used  to  predict  bubble  distributions;  its  primary  utility  is  as  an  analytical 
tool  for  investigating  relationships  among  the  various  complex  factors  that 
affect  the  distribution  and  for  refining  models  of  the  upper  ocean.  Future 
work  will  be  concerned  with  analysis  of  more  complex  ocean  models. 
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APPENDIX 

It  is  shown  here  that  the  Jacobian, 

Y(tft.,.ag£|eiL  ,  (A_1} 

3(r;v;£') 

used  to  transform  the  volume  element  d7x  between  two  points  on  a  characteristic 
curve, 

d3rd3vd£  =  YCtft'JdVdVdA',  (A-2) 

satisfies  the  differential  equation 

dY/dt  =  Zty  .  (A-3) 

Using  the  boundary  condition  y(tjt')  =  1,  this  equation  can  be  integrated  to 
obtain  t 

Y(t,t')  =  exp[7  y 


(t,f)  =  exp[7  Zt^'VVV^'l 


where  the  indicated  integration  is  along  the  characteristic  curve  passing  through 
(r;vj£')  and  (r,v,Jl) . 

Referring  to  the  notation  (5)  ,  define  the  elements  of  a  7  x  7  matrix 

Yij  =  3Si/a9j''  i'J  =  1/2,..., 7  (A-5) 

then  the  Jacobian  is  the  determinant  of  this  matrix, 

Y  =  det  ( (Yij ) )  =  J^Z    (-D 1+j Yi jci  j  >  (A"6) 

(or  j) 

where  C.  .  is  the  "minor"  of  the  element  y.  .  (C  .  is  y  with  the  i   row  and 
i:  ID   ^0 

3th  column  deleted) . 
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If  the  following  elements  are  defined, 


*!  ~=   V  R2  E  V  *3  '=   V  R4  E  ax 


(A-7) 


then  the  characteristic  equations  (11)  can  be  written 

dq±/dt  =  Ri(q1/q2/...,q7/t)   i  =  1,2,... ,1  (A-8) 

Taking  the  partial  derivatives  of  these  equations  with  respect  to  the  q . '  and 
using  the  chain  rule,  the  set  of  differential  equations  for  the  y.  .  is  obtained: 

dYiJ  =  9Ri  S^^L    Yw  *•'$'   "  1'2 7        <A"9> 

3F~   ^qT   Z-,9qv    D 
D     k   K 

These  equations  are  used  to  expand  the  derivative  of  (A-6)  taken  with  respect 
to  the  parameter  t: 


i+j  dy.  •      ^"~ *    i+D  "T""1  3R. 


dy  ^"^  ^"^    i+j  dy....      ^— '    i+j  ^ — '  3R 

dt  Z. >  £ — '  13  c — t  *qu 

ID  i,d         k 


Using  (A-6) ,  this  expression  can  be  rearranged  to  yield  the  following: 


r — •     i+j     f  3R.      C"^  3R.     -j 


dy   ^ — •    i+j 

at 

i/j  A      k^i 


(A-ll) 


^— »3R.   r — *     i  8R.  r-* 


l)  y,  •  c.  . 
k]  id 
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The  term  2. ►  (~1)  Yv-  c-  •  ^  zero  since  k^i.  This  can  most 

easily  be  seen  by  expanding  each  determinant  C.  .  by  minors  of  the  elements 

Ykm'  "*> 

Since 


L 


9qi 

1 


3R.  .     — 

_i  =  ^   .  t  +  |v  -  Z  (A-12) 

3q.     v       3£      t 


the  result  (A-2)  is  obtained. 
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LIST  OF  FIGUEE  CAPTIONS 

Figure  1: 

The  velocity,  v,  of  bubbles  starting  at  rest  and  rising  from  an  initial 

depth  z  =  -20  meters.  The  initial  radius  is  I   ,  (in  microns) ,  the  terminal 
v        o  o 

velocity  is  vT,  and  (z-z  )  is  the  height  above  the  release  point. 

Figures  2,  3 ,  and  4  : 

Bubble  radius  in  microns  as  a  function  of  depth  in  meters  for  various 
initial  radii,  I     (at  z  =  -20  meters) , calculated  using  three  different  expres- 
sions for  the  terminal  velocity,  v_.  The  solid  line  is  I   ,  obtained  by  inte- 
grating (38)  with  v  =  v_  (equation  42),  the  broken  line  is  H—  (equation  43a),  and 
the  dashed  line  is  £_p  (equation  44a) .  The  gas  diffusion  model  used  (Wyman,  et  al 
1952)  assumes  a  uniform  partial  pressure  of  one  atmosphere,  f  =  1,  and 
RTGK/4tt  =  10" 6  meters/second. 
Figure  5  : 

The  distribution,  $  (z  ,  £  ),  in  relative  units,  of  bubbles  rising  off  the 
bottom  of  Saanich  Inlet,  B.  C.  (McCartney  and  Bary  1965)  where  z  =  -197  meters 
and  the  bubble  radius,  I   ,  is  in  microns. 

Figures  6,7,8,  and  9: 

The  distributions  $(z,£)    (solid  line)   and  $~(z,£)    (dashed  line)   in  relative 
units  at  the   depths    -172m,    -147m,    -122m,    and   -97m  where  bubble      radius, 
i,    is    in  microns.       $      neglects   gas   diffusion  while    $    assumes   a  uniform 
partial  pressure   of  one   atmosphere,    f  =    1,    and  RTGK/4tt  =   10  meters/ 

sec    (Wyman,    et    al    1952) . 
Figure   10: 

The  distribution,   $-(z,£),  in  relative  units,  vs  radius,   I,  in  microns,  for 
z  =-47m     and  z  =   o. 
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Figures  11  and  12: 

The  distribution,  $(z,£)  (solid  lines),  in  relative  units,  vs  radius,  I, 
in  microns,  for  various  depths,  in  meters,  from  a  distributed  source,  s(z,£) 
(dashed  lines;  equations  61  and  62)  .  The  same  gas  diffusion  model  used  to 
calculate  $  in  Figure  5  is  used  in  this  example. 
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